A perturbative approach to the Bak and Sneppen Model 



X 



G. Caldarelli\ M. Felici^'^, A. Gabriclli^ and L. Pictronero^ 
^ INFM - Umta di Roma 1 "La Sapienza", P.le A. Mora 2, 00185 - Roma, Italy 
^ Lah. PMC, Ecole Polytechnique 91128 Palaiseau, France 

We study here the Bak and Sneppen model, a prototype model for the study of Self-Organized 
Criticality. In this model several species interact and undergo extinction with a power law distribu- 
tion of activity bursts. Species are defined through their "fitness' whose distribution in the system 
is uniform above a certain universal threshold. Run time statistics is introduced for the analysis 
of the dynamics in order to explain the peculiar properties of the model. This approach based on 
conditional probability theory, takes into account the correlations due to memory effects. In this 
way, we may compute analytically the value of the fitness threshold with the desired precision. This 
represents a substantial improvement with respect to the traditional mean field approach. 
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I. INTRODUCTION 

In the recent years, many models mimicking the scale 
free behavior exhibited by natural phenomena like river 
basins , fracture dynamics |^,|| , earthquakes j|] , have 
been extensively studied. The main features of these 
models are the lack of spatial and temporal character- 
istic scales, the evolution through intermittent bursts of 
activity and the absence of the fine-tuning of some pa- 
rameter to reach the critical state. To describe all this, 
the concept of Self-Organized Criticality (SOC) |^ has 
been introduced, and many works have been devoted to 
clarify its real nature . 

We focus here on the evolution model introduced by 
Bak and Sneppen (BS) in 1993 §. This model is the 
prototype of a wide class of SOC models, characterized 
by a deterministic dynamics in a medium with quenched 
disorder. The quenched noise is independent on time and 
represents the disordered environment where the system 
evolves. BS model is defined by a discrete set of L species 
arranged on a one dimensional network. Each species i 
is defined by a "fitness" Xi given by a real number in the 
interval (0,1). Time evolution is discretized and at any 
time step the species with the lowest fitness is removed 
from the set. At the same time also the two species neigh- 
bors are removed. This should model a food web where 
extinction of one species afi'ects also the survival of preda- 
tion and prey species. Three new species with randomly 
extracted new fitness enter then the system. Hereafter, 
we shall call this process an "updating" of the species 
fitnesses. After a transient period, the system reaches 
spontaneously a stationary state characterized by two 
main features: (i) the fitnesses are uniformly distributed 
between a threshold value Xc and 1; (ii) the dynamics 
evolves as a sequence of critical avalanches |^ , whose du- 
ration s is power-law distributed [?]: P{s) ^ s^"^, where 
T = 1.07 0. 

To study in an analytical way this kind of processes, 
a method based on conditional probability has been re- 
cently introduced. This method, called Run Time Statis- 
tics (RTS), provides a powerful tool to study how the 
system stores information on the disorder during its evo- 



lution. This method has been also applied to the class 
of models derived by the Invasion Percolation [ p^|Jl^ 
to compute the asymptotic behaviour of the histogram 
equation. We apply here RTS to the BS model in order 
to find an analytical approach to solve the model. 

As regards the analytical results available for the BS, 
the only successful approach has been the mean field ap- 
proximation. Scale relations and an equation describ- 
ing the hierarchy of the avalanches ||l^ , allow to reduce 
the number of independent exponents to one. A kind 
of e-expansion was also introduced to compute the crit- 
ical exponents by performing an expansion around their 
mean field value. This novel approach allows to compute 
in a perturbative way the value of the critical threshold 
linking it to the avalanche exponenti r. In the following 
sections we are going to introduce the main features of 
the model, the basic computation of the run time statis- 
tics and the result of this approach with respect to the 
BS. A preliminary letter with some of the results has al- 
ready been published in Ref. p5[ |, here instead we are 
going to fully develop the derivation of such results. 



II. CRITICAL DYNAMICS 

The dynamics of the BS model can be viewed as a 
branching process of branching ratio L, if L is the num- 
ber of species in the system. At every time step, one of 
the L quenched numbers assigned to the sites is selected 
to be updated. The same is done for the two nearest 
neighbors. It is then possible to represent this process 
as a tree-like picture where every node represents a state 
from which the system can reach L possible states at 
the subsequent time step. Therefore at time t there are 
L* possible states that are reached through different 
paths. Because of the deterministic nature of the model, 
time evolution is determined by the realization of the dis- 
order {x}. This means that if we know the initial set of 
values Xi, i — 1, . . . , L, and the three numbers extracted 
at each time-step then we also know a priori the evolu- 
tion of the system. Otherwise, to give a description of 
the behavior of the generic system, we should consider 
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a statistical average over the possible realizations of the 
quenched disorder. RTS provides an iterative algorithm 
to assign to each path (sequence of events) its statistical 
weight according to the laws of conditional probability. 

We consider the system in the critical state. In this 
stationary critical state of BS model almost all the fit- 
nesses lie above a threshold Xc- The distribution function 
of the fitnesses in the critical state, called the histogram 
is found to be uniform above Xc, and almost zero 
below Xc- The critical state is characterized by a power 
law distribution of the duration s of critical avalanches: 
P{s) ^ s^'^. Critical avalanches are defined as a sequence 
of events Xmin[s) < Xc, where Xmin{s) is signal, i.e. the 
value of the minimal fitness at time s [D . Then a criti- 
cal avalanche begins each time the signal reaches a value 
larger than Xc- Because of the shape of the histogram, 
that is of order for x < Xc and constant for x > Xc- In 
the critical state the signal reaches at most values near to 
Xc and the dynamics is a sequence of critical avalanches 
[^. Then the site giving rise to a new avalanche, called 
the initiator, has a fitness near to Xc.This means that 
when the initiator is selected, all the other sites have fit- 
nesses larger than Xc- 

Since the critical avalanches are independent each 
other 1^, we can consider the dynamics within one 
generic critical avalanche. This is then a fair descrip- 
tion of the dynamics of the system in the critical state. 
Let us define the set of active sites At as the sites covered 
by the avalanche, i.e. the sites whose fitness has been up- 
dated at least once since the beginning of the avalanche. 
If we consider the dynamics within an avalanche, the pos- 
sible events at time t regard only the sites in At because 
the selection of any other site would imply the end of 
the avalanche. Indeed, sites not belonging to At have a 
fitness larger than Xc (by definition of At their fitness re- 
mained unchanged since the first step of the avalanche). 
Therefore the evolution of the avalanche can be seen as a 
branching process where the branching ratio is not fixed. 

Without loss of generality, we take the time origin at 
the beginning of a critical avalanche, and the origin of the 
coordinates in the initiator site. At t — the initiator is 
the site with the smallest fitness (i.e. extremal rule) and 
its quenched variable (in the stationary state) has a value 
near to Xc- Since we are assuming that the stationary 
state in the system is the critical one, all the quenched 
numbers in the system are distributed following the sta- 
tionary distribution ^(x), apart from corrections of order 
j^. Updating at time t = Q affects the initiator and the 
two nearest neighbors. At—i is then composed by the 
three sites { — 1,0,1}. The three variables corresponding 
to these sites are distributed following the uniform prob- 
ability density fo{x). If the avalanche proceeds, there are 
three possible events leading to three different configura- 
tions for the system (see Figjl]) . Then a path of length 
2 can be realized in the following three ways: i) growth 
of the initiator at time t = and then growth of its left 
neighbor at time t = 1; ii) growth of the initiator at 
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FIG. 1. Diagrammatic plot of the first three steps in an 
avalanche tree. The initiator i = is selected at t = 0. 
At each step one has to consider all the possible offsprings. 
Non-active sites are represented by a filled circle; updated 
active sites are represented by with a simple cross. Crossed 
empty circles represents instead active sites not updated at 
the previous time steps. 
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TABLE I. Number of possible paths as a function of time. 



time t = 0, and successive growth of the initiator at time 
t = 2; iii) growth of the initiator at time t = and then 
growth of the right neighbor at time t — 1. As the length 
t (i.e. the number of steps) increases, the number Nt of 
possible paths of that length increases fastly. For exam- 
ple, avalanches lasting two time steps can occur in the 
previous three ways, but there are eleven ways to form 
an avalanche of three time steps (see Fig.^. In table || we 
report the number Nt of possible paths whose length is t. 
It can be shown that Nt grows roughly as tl. Given this 
picture of the stationary state, it is evident that a de- 
scription of the model evolution can be achieved through 
a description of the growth paths. Run Time Statistics 
(whose approach is described in the next section) helps in 
sorting out the most probable paths, thereby extracting 
the available information on the process. 



III. THE RTS APPROACH 

Through the RTS we are able to compute the statis- 
tical weight of all the possible paths corresponding to a 
critical avalanche of a fixed duration 17|. The RTS 
has been developed in order to extract the maximal in- 
formation available from the knowledge of the dynamical 
history followed by the process. The information is stored 
in effective probability density functions (PDF's) for the 
variables {x}. These effective density functions at a cer- 
tain time t, are used to compute the conditional proba- 
bilities of all the possible events at time t + 1 given the 
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state at time t. The probability of a sequence of events 
(that is a fixed path), is then factorized in the product of 
these one-step-probabihties. The time dependent PDF's 
are obtained by applying the laws of conditional proba- 
bility. At the beginning the only information available on 
the disorder is the probability density fo{x) from which 
the quenched variables are extracted. The conditional 
probability laws in the following steps modify the shape 
of the probability distribution once the previous history 
is known. 

In the case of BS, at every time step the smallest num- 
ber is removed from the system. It is then intuitive that 
if the lifetime of a species, i.e. the number of tests the 
species survived in the search of the minimum, is large, 
its fitness is (probably) also large. Henceforth, the longer 
the lifetime the higher the probability to deal with a large 
value of the fitness x, and the smaller the probability to 
be removed. This concept can be easily formalized by 
using the theory of conditional probability. 

Let us introduce the time dependent PDF's fi,t{x) giv- 
ing the probability density of the quenched variable of 
site i at time t. If mi^t{x) is the probability density func- 
tion of fitness at time t (assuming it corresponds to site 
i), we have that mi^t{x)dx is the probability that the fit- 
ness of site i has a value between x and x + dx (given 
that at time t of a fixed path, i is the minimum with fit- 
ness smaller than Xc] this last condition ensures that the 
system is still under the same critical avalanche). Then 
nii^tix) is given by: 

1Tli,t{x)=P{x < Xi < X + dx\Xi — min{x^ At-n < ^c) = 

P [x < Xi < X + dx n Xi — min{x}At-i, < ^c) 
P (xj = min{x}At-i, < ^c) 

J[^::^f^,t~l{x)ll^^^^_^_^^y dxkfk,t-iixk)x < Xc 

X> Xc (1) 

Where we have defined the one step probability Hi^t- 

Hi,t = P{xi = min{x}At,Xt < Xc) = 



This is because of the condition that the minimum fit- 
ness is less than Xc- The complementary probability 
(l — X^ieAt /^«,*) ^'^ then the probability that the mini- 
mum is larger than x^, that is that the avalanche stops 
at time t. 

In a similar way we can be obtain the probability den- 
sities fk,t{x), with k G At-i'. 

fk,t{.x)^P[x < Xk < X + dx\xi = min{x}At-i , Xi < Xc) = 



P [x < Xk < X + dx n Xi = min{x}At^i, < Xc) 
P (xi = min{x}At-i,Xi < Xc) 

j^fkA^) lo hAxi)dxr Ylj dxjfjAxj)x < Xc 

■j^JkA^c) Jo" fiAxi)dxi Ylj dxjfjAxj)x > Xc 

(4) 

where j e At-i — {i, k}. In this way the effective PDF's 
are conditioned to the whole history from time (be- 
ginning of the avalanche) to time t because of the step- 
by-step algorithm through which they are obtained. We 
also notice that if the minimum fitness is less than Xc, 
this implies that /i,t(x) = fiA^c) for x > Xc- As it has 
been pointed out in [p|, the dynamics involves only the 
quenched numbers below Xc, this is the reason why the 
system does not acquire information on the variables in 
the region x > Xc- 

These formulae for the effective probability densities 
hold if we assume that the probability density of the 
whole set of variables Ft{{x}) could be factorized in the 
product of the one-variable probability densities at any 
time. Actually this is true only at time t ~ 0, when the 
{x} are uncorrelated, while at later times the extremal 
dynamics induces correlations among them. Nevertheless 
the approximation 



dxfi^t-i{x) Y]_ / d^jfj-t'iixj)- (2) 



(li^t represents the probability that site i has the minimal 
fitness at time t (smaller than Xc), given the path followed 
up to time t — 1 (all the information about the past 
steps is included in the effective PDF's fjAx))- In both 
eq.(||) and eq.(||) we consider only the fitnesses {a;} At of 
the active sites At because the others actually do not 
participate in the dynamics, being larger than Xc- It is 
important to notice that 



Ft{{x})^l[f,Ax) 



(5) 



has proven to lead to results in good agreement with data 
from simulations |jl^, and to give rise to an error which 
is negligible for large values of the system size L . 

Once the one-step-probabilities /i^.f have been com- 
puted, the probability of a given path Ct (i.e. a fixed 
sequence of events from time to time t) is given by: 



(6) 



t'=i 



(3) 



This probability is the probability of the path Ct aver- 
aged on the disorder. 

Let us consider what would be the rigorous computa- 
tion of this quantity. The extremal rule can be formu- 
lated by defining a growth probability r]i{{x}) for site i 
to be selected, given by: 
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FIG. 2. Statistical weights ^(Ce) of the paths Ce of length FIG. 3. Effective probability density of the site 1 at differ- 
6. Paths are numbered from 1 to Nq following the order in a ent times of the path Ia2a3a4a. The probability densities are 
tree-hke diagram like the one shown in Fig|l[ labeled with the indexes r, t giving the "age" of the quenched 

number, and the time step of the path. 



1 Xi = min{x} 
Xi min{x} 



(7) 



Then the probability of a given path Ct is given by: 



(8) 



which can assume only the values 1 or 0. The exact com- 
putation of the mean probabihty of the path Ct consists 
in taking the average over the reahzations of the disorder. 
This average cannot be factorized: this means that it is 
not possible to obtain the mean probability of a path Ct 
by simply multiplying the mean probability of the path 
Ct-i, corresponding to the first t — 1 steps of path Ct, for 
the mean probability of the last step. On the contrary, 
the weight W{Ct) is factorized as in eq.(^) because the 
one-step probabilities fii^t are conditional probabilities. 
The weights W{Ct) of the paths of length 6 are plotted 
in Fig.(||). In Fig(g) are plotted the effective probability 
densities fi.t{x) of a site i for different values of t: the 
site chosen is not updated from time t — I to time i = 4, 
the corresponding PDF is then modified in such a way 
to approach the histogram shape. 



IV. THE HISTOGRAM EQUATION 

We now introduce an equation for the histogram ^t{x), 
that represents the average probability density of the 
quenched variables at time t. In the limit i ^ oo we 
obtain an equation for the stationary histogram. Let us 
introduce the function ht (x) defined in the following way: 



ht{x) = L'S>t{x) 



(9) 



where L is the size of the system. Then ht{x)dx is the 
average number of quenched variables in the system in 
the interval [x, x + dx\ at time t. Since the species up- 
dated at each time step are the minimal one with the two 
nearest neighbors, we can write a balance equation : 

ht+i{x) = ht[x) - mt+i[x) - {fi{x) -t- f2{x)) + 3/o(a;) 

(10) 

where mt{x) is the probability density of the minimal 
variable; fi,2{x) are the probability densities that the 
two nearest neighbors variables would have at time t+l 
if they were not updated. The FDF's fi,2{x) are given 
by the first line of eq.(^ with k ~ i — I and k = i + I 
respectively. Because of the self averaging property of 
this function the result obtained by taking the asymp- 
totic limit coincides with the one obtained by averaging 
over the possible realizations of the disorder: 



lim ht{x) = L^{x) = L{^{x))m 



We then obtain: 



(mix)) + ihix) + f2ix)) 







(11) 



(12) 



where fo{x) = 1. 

To compute the average {fi{x) + f2{x)) we use the 
weights W{Ct) obtained by applying the RTS algorithm: 
(/i(x) -I- f2{x)) is given by averaging over the paths Ct, 
weighed with the W{Ct) 



if lix) + f2{x)) 



EZiEcWiCt) f?\x) + f^\x) 



fCt, 



EZiEc. w{Ct) 



(13) 
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To compute this quantity to the order n, we perform an 
exact enumeration of the paths of length t < n. For each 
path Ct, f?2{x) are computed by iterating the formu- 
lae in eq.(Q) to obtain the effective probability densities 
at time i + 1 of the nearest neighbors sites of the site 
selected to grow at time t of the path Ct- The sum in 
eq.(|l^) contains also terms proportional to <^{x). This 
happens when one of the nearest neighbors of the ex- 
tremal site does not belong to the set of active sites, 
its probability density is (^{x). Moreover Xc appears ex- 
plicitly in eq. (p^) because both the probability densities 
fi,2{x) and the weights W{Ct) depend on it, (since Xc 
has been introduced as a parameter by imposing that 
the system is under an Xc-avalanche). 

The minimum probability density is averaged in a dif- 
ferent way, following a mean field argument. We now 
consider the Generalized Run Time Statistics (l9| , which 
gives the correct form for the effective probability densi- 
ties in the case of a stochastic dynamics. Let us suppose 
to have a growth probability r]{{x}) depending on the 
quenched disorder different from extremal rule. Also in 
this case the system stores information during the evo- 
lution. The eq.^ must now be modified in order to take 
into account the probability r]i{{x}). The one-step prob- 
ability is obtained in the following way 



/ fi,tixz)dxt ■■■ Yi ^^j/j-t(^i)^«({^}'4t) 

Jo Jo Jo . r,-i 



(14) 



where now At is given by all the sites in the system. The 
probability density of the minimum site is consequently 
given by: 

1 

m^,t+\{x) ^ / dxib{xi -x)fi^t{xi)- 

f^i,t Jo 



■ ■■■ T\'ni{{x}At)fk.t{xk)dxk (15) 
Jo Jo J, 

where it is the quenched number of the site selected at 
time t. 

We now consider the growth probability given by: 



e T 



(16) 



e T 



In this way in the limit T — )■ we recover the extremal 
rule: 



lim 



e T \ I Xi = 'min{x} 

^ Q-^- I Xi min{x} 



(17) 



We then substitute the expression eq.(p^ in eq. dH) to 
realize the average over the disorder and then the limit 
T ^ 0. 



We average on the paths Ct of length t and then con- 
sider the limit t ^ oo. The average on the paths Ct is 
taken by averaging first on the last step and then on the 
paths Ct„i: 



{m^^t{x))ct = (("^^.t(a;))^,)ct-l 
The first average gives: 



(18) 



{m^^t{x))^^_, = y2^^''*{ —fi:t{^)- 



/"•••/" \\dxjfj,t[x)—^ 

Jo Jo , e T 



e T 



Then, averaging over the paths Ct — 1 we obtain: 

((m,,t(a;)),Jct-i = 



Jo 



e T 



e ^ + E 



e T 



Jo 



Yldxjifj^tix)) — 



e T 



(19) 



Taking the limit t — > cxo we have (fj.t{x)) = $(a;): 



(m(a;)) = L$(a;) 



Jo 



Y[dxj^ix:j) 



e T 



(20) 



Let us consider the integral in the RHS of eq.(pO|): it is 
the average of the function g{{xk}) 



which is a function of the (L — 1) variables {xk}k^i ■ 

If we put Zj = e~~ and Z = X^fc^fc; according to the 
central limit theorem the deviation from the mean value 
{Z) = {L — ^){z) is negligible. We can then make the 
approximation {F{Z)) « F{{Z)): 



Jo 



X{dx,^{xj) — 



e T 



e~T + (L - l)(e"T") 



We now introduce the parameter Xc by defining (z) 
(e~T) = e~T*. We finally obtain: 



{m{x)) = L^{x) 



1 + {L- l)e ''^■r 



(21) 
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We can see that if we assume for ^{x) the following be- 
havior: 



$(x) 




(22) 



then the average of the minimum density function given 
above has the expected behavior, in the limit T —^ 0: 



(m(x)} 



0(1) 




(23) 



Thus the parameter Xc introduced above is the critical 
threshold in the histogram function. 

We can now turn to the eq.(^2|). After use of eq.(|21 
and eq.(0|), the equation takes the following form: 



1 



— (Xf. — x) 

1 + {L-I)e^^ 



A^"'>{xc)^{x) +B^"'>{x,Xc) 
Di^^xc) 



-3 = 



(24) 



where the coefficients A'-'^^Xc), i3'"''(a;, Xc ) and ^("^(xc) 
are given by eq.(|l3|) with truncation at order n. There is 
only one value of Xc for which ^{x) is normalized. This 
value gives the value of the BS threshold. 

It is easy to verify that the first order of approxima- 
tion, n — 1, corresponds to the mean field approximation. 
Indeed, the first order of eq.(n3h is: 



(25) 



{fi{x)+h{x))^^^ ^2<S>{x) 



This is essentially the random neighbors assumption 
since in this way all the correlations among the species 
in the avalanche are neglected. By substituting eq.(|2F 
in eq.(p3) we obtain: 



$(a;) = 



(26) 



l + (L-l)e 

that in the limit T ^ 0, becomes: 



$(a;) 



Then, imposing the normalization condition 



^{x)dx = 1 



(27) 



(28) 



we obtain |(1 



- Xc) — 1, verified for 



value obtained in the mean field case ll^ 



■i , which is the 
. It is worth to 
notice that in the eq.(27) we have obtained analytically 
the behavior <I>(x) ^ j- under Xc- 

We solved numerically eq.(|lj) from order n = 2 to 
order n = 7. The values of Xc obtained for different values 
of n are plotted in Fig.0. By considering part of the 



correlations among species Xc becomes larger than the 
one obtained in the mean field approximation ( a;c = 1/3 
). The best evaluation is Xc{n — 7) ~ 0.465, larger than 
the mean field result but still quite far from the value 
obtained from simulations. Nevertheless, it is possible 
to verify that the behavior of Xc{n) is compatible with 
an asymptotic value Xc{n oo) ~ 0.66. We made a 



fit with the fitting function xdn) = 0.66 



which is 



found to be well compatible with the given asymptotic 
value (see Fig.|). The fit values are a = 0.291 ± 0.003 
and b = 0.20 ± 0.03. The small value of b is due to 
the fact that the avalanche duration distribution P{s) is 
characterized by a small exponent (r ~ 1-07), henceforth 
all the sizes s are important for statistics. 

One can use Xc{n — > oo) to evaluate both the 
avalanche exponent r and the average minimum distri- 
bution m^x) = {mi{x)). The exponent r can be obtained 
from the function VL{t) = X^Ct ^i^t) for t ranging from 1 
to the maximal possible n. This function is proportional 
to the probability that the avalanche lasts at least t time- 
steps. Thus, in the scaling regime Vl{t) ~ t~'^^^. Making 
this hypothesis, and substituting the value Xc{n —>■ oo) 
in the expressions giving the weights W{Ct), one finds 
T{n — 7) ~ 1.05 (see Fig.^), which is in agreement with 
the known numerical value. 

Finally, we can obtain an approximation of the proba- 
bility density function of the minimal fitness m{x) in the 
stationary state. Let us consider the eqs.^ and If 
we take the limit T ^ in eq.01^ we obtain: 



{■m{x)) = L^{x)6{xc ~ x) 



(29) 



If we now turn to eq.(p4[) and we solve it for $(a;) us- 
ing the value Xc{n — > oo), we obtain the function m{x) 
reported in Fig.^. In the same figure, this result is com- 
pared with the numerical distribution of the minimal fit- 
nesses obtained in the numerical simulations. In spite of 
the strong approximation (the paths considered are only 
those of length < 7), the agreement is quite good. 

In conclusion, this paper presents a perturbative ap- 
proach to the BS model, based on the probabilistic frame- 
work called Run Time Statistics. The detailed derivation 
of the self-organized threshold Xc, the avalanche exponent 
T, and the stationary distribution of minimal fitnesses 
m(x) is presented here. Through RTS we are able to 
improve the agreement between the numerical and the 
theoretical values found for this model. 



[1] I. Rodriguez-Iturbe and A. Rinaldo Fractal River Basins, 
chance and Self- Organization, Cambridge University 
Press, Cambridge UK. 

[2] G. Caldarelli, A. Giacometti, A. Maritan, I. Rodriguez- 
Iturbe and A. Rinaldo Phys. Rev. E. 55, R4865 (1997) 



6 



0.48 




n 



FIG. 4. Empty points represent the values of Xc obtained 
by the application of the RTS algorithm from n = 2 to n = 7. 
The continuous lino represents the fit curve Xc{ri) = 0.66 — aa::'' 
with a = 0.291 ± 0.003 and b = 0.20 ± 0.03. The insert shows 
the behavior of f2(n) up to n = 7. Assuming f2(n) ~ n"''"'"^ 
as a good approximation also for small n, one finds t ~ 1.05. 
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plication of the RTS algorithm, considering all the possible 
avalanche paths up to a maximal time n = 7 after the selec- 
tion of the initiator, and assuming Xc{n —^ oo) = 0.66. The 
points gives the numerical behavior evaluated in extensive 
simulations. 
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